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Within a dipolar Poisson-Boltzmann theory including electrostatic correlations, we consider the 
effect of explicit solvent structure on solvent and ion partition confined to charged nanopores. We 
develop a relaxation scheme for the solution of this highly non-linear integro-differential equation 
for the electrostatic potential. The scheme is an extension of the approach previously introduced 
for simple planes (S. Buyukdagli and Raff Blossey, J. Chem. Phys. 140 , 234903 (2014)) to nanoslit 
geometry. We show that the reduced dielectric response of solvent molecules at the membrane walls 
gives rise to an electric field significantly stronger than the field of the classical Poisson-Boltzmann 
equation. This peculiarity associated with non-local electrostatic interactions results in turn in 
an interfacial counterion adsorption layer absent in continuum theories. The observation of this 
enhanced counterion affinity in the very close vicinity of the interface may have important impacts 
on nanofludic transport through charged nanopores. Our results indicate the quantitative inaccuracy 
of solvent implicit nanofiltration theories in predicting the ionic selectivity of membrane nanopores. 
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I. INTRODUCTION 

Electrostatic interactions in confined liquids are rele¬ 
vant to various nanoscale processes ranging from macro- 
molecular interactions ffi m to nanofiltration 012 and 
nanofluidic transport 01- An accurate theoretical for¬ 
mulation of these interactions has been a major chal¬ 
lenge in the field of theoretical soft matter physics. The 
description of electrostatics in nanoscale systems has 
been mainly limited to Poisson-Boltzmann(PB) formal¬ 
ism 0 0 known to suffer from two major limitations. 
First, the PB equation provides a mean-field (MF) de¬ 
scription neglecting electrostatic correlations resulting 
from many-body interactions between charges. The sec¬ 
ond weak point of the PB approach is the dielectric 
continuum approximation bypassing the extended charge 
structure of solvent molecules hydrating the ions of the 
solution. 

In confined systems where the distance between 
charged macromolecules approach the Bjerrum length 
i w ss 7 A, the above-mentioned limitations become dras¬ 
tic. First of all, the dielectric contrast between the 
macromolecule of low dielectric permittivity and the elec¬ 
trolyte of large permittivity is known to induce surface 
polarization effects. These non-MF effects that can¬ 
not be captured by the PB formalism are at the origin 
of i) attractive van-der-Waals forces stabilizing macro¬ 
molecules [T] and ii) dielectric exclusion phenomenon set¬ 
tling the salt selectivity of membrane nanopores in cells 
as well as in artificial nanofiltration devices 0. Then, 
by ignoring the interaction between solvent molecules 
and the membrane, the PB formalism assumes the same 
solvent density and permittivity in the pore and the 
bulk reservoir. This dielectric continuum approximation 
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is expected to result in two artefacts. First, solvent- 
membrane interactions are known to induce repulsive 
image-solvent effects leading to a reduced pore dielec¬ 
tric permittivity [9]. This implies that solvent implicit 
electrostatic theories [Tf)l ITT] overestimate the pore per¬ 
mittivity and the strength of image-charge forces exclud¬ 
ing ions from the nanopore. However, the same dielectric 
contrast between the pore and the reservoir should lead 
to an ionic Born energy barrier, i.e. an additional force 
favouring ionic rejection from the pore medium. Hence, 
the dielectric continuum limit giving rise to opposing 
artefacts is an uncontrolled approximation that has to 
be relaxed by the explicit consideration of the solvent 
charge structure. 

Improvements over the mean-field and dielectric con¬ 
tinuum approximations have been mainly considered sep¬ 
arately. Within the dielectric continuum limit, elec¬ 
trostatic many-body effects have been included at one- 
loop [T2UT6] and variational levels [TUI HU fTTHST] . In 
the MF approximation, explicit solvent theories have 
been developed by modeling water molecules as point 
dipoles [221 S3]. The point dipole approximation be¬ 
ing unable to account for non-local electrostatic inter¬ 
actions, we introduced in Ref. [22 the first solvent ex¬ 
plicit formulation of non-local electrostatics with finite 
size dipoles. We explored next the MF non-linear re¬ 
sponse regime of this model in Ref. [25]. At this point, 
one should also mention the phenomenological non-local 
electrostatic theories developed in Refs. [55H50] . Within 
the point dipole approximation, electrostatic correlations 
were considered in bulk liquids in Ref. [3T] and at single 
charged planes in Ref. [52] in order to scrutinize the ef¬ 
fect of interfacial solvent configuration on the differential 
capacitance of low dielectric materials. 

The unification of non-locality and electrostatic corre¬ 
lations was introduced first in Ref. [9]. In this previous 
work, we derived the solvent-explicit variational equa¬ 
tions with finite size dipoles. We solved these equations 
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for an ion free solvent confined to neutral slits and also 
for an electrolyte in contact with a single charged inter¬ 
face. In the present article, we take a step further and 
develop a numerical scheme for the solution of the non¬ 
local PB equation (Eq. Q in the main text) in charged 
membranes of slit geometry. In other words, the consid¬ 
eration of the liquid confinement and membrane charge 
is the novelty that had not been treated together in our 
previous article [5]. Within the framework of this for¬ 
malism, we investigate the effect of surface charge and 
pore size on the configuration of solvent molecules, the 
liquid charge, and ion densities in the nanoslit. In par¬ 
ticular, we show that due to an amplified surface field 
induced by the reduced solvent dielectric response at the 
interface, an ionic structure formation characterized by a 
double counterion concentration peak takes place. This 
structure formation corresponding to an enhanced coun¬ 
terion affinity is the most relevant result of the present 
work. We emphasize that the formalism developed herein 
provides the first explicit solvent theory of membrane se¬ 
lectivity. The limitation of our theory and possible im¬ 
provements are discussed in the Conclusion. 



FIG. 1: (Color online) Schematic representation of the 
nanoslit with thickness d and negative surface charge — a s . 
Dipoles of size a (red) hydrate cations (yellow) and anions 
(blue). The dipolar orientation with respect to the a—axis is 
characterized by the projection a z . 


II. NON-MF NLPB FORMALISM 

In this part, we review briefly the explicit solvent the¬ 
ory of electrolytes previously developed in Ref. |5j. The 
charged liquid is composed of solvent molecules and point 
ions confined to a slit pore. The composition of solvent 
molecules and the pore is depicted in Fig. |T| The pore is 
composed of two infinitely large and negatively charged 
planes located at z = 0 and z = d. The charged particles 
in the pore are in chemical equilibrium with an external 
charge reservoir. Each solvent molecule is modeled as 
a finite-size dipole composed of two oppositely charged 
point particles ±Q separated by the fixed distance a. The 
reservoir concentration of solvent molecules is p s b . There 
are p ionic species in the solution and each species i has 
bulk concentration pn, with i = 1, ...,p. 

The beyond-MF NLPB equation for the electrostatic 
potential <fi(z) in the solvent-explicit liquid was derived 
in Ref. [9]. This equation of state reads 


k T 

-d z £ 0 {z)d z (/)(z ) + Y QiPi(z) + Q [Ps+{z ) - Ps— (~)] 
= ( 1 ) 

where ks is the Boltzmann factor, T = 300 K the am¬ 
bient temperature, and e the electron charge. The per¬ 
mittivity function is given by the piecewise form £q(z) = 
Eo9(z)9(d — z) + £ m [0(— z) + 9(z — d)], with £ 0 the vac¬ 
uum permittivity and e m the membrane permittivity. 
In the present work, dielectric permittivities will be ex¬ 
pressed in units of the vacuum permittivity. Moreover, (/,; 
stands for the valency of the ionic species i, and the sur¬ 
face charge density reads a(z) = — a s [£( 2 ) +6(d— z)\. 
Finally, the ion number density and the density of the 


positive and negative charges on the solvent molecules 
are respectively given by 


Pi(z) 

Ps±(z) 
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In Eq. the integration variable a z corresponds to the 
projection of the dipolar alignment vector a onto the z 
axis (see Fig. T]) and the integration boundaries taking 
into account the rigidity of the interfaces are given by 


a\(z) = — min(a, *) (4) 

a 2 (z) = min (a,d—z). (5) 

We introduced in Eqs. ([ 2 ]) and ([3| the renormalized 
ionic and dipolar self-energies Svi(z) and 5vd(z, a z ) whose 
explicit form is reported in Appendix [A] The relaxation 
algorithm that solves the integro-differential equation (|T|) 
is explained in detail in Appendix [b] ^33). We note that 
this algorithm generalizes to confined slits the scheme in¬ 
troduced in Ref. |9] for the single interface geometry. Fur¬ 
thermore, we note that in the present work, the charges 
composing the dipolar solvent molecules will be monova¬ 
lent, that is Q = 1. The solvent molecular size will be set 
to a = 1 A, and the bulk solvent density will be taken as 
the liquid water density p s b = 55 M. Finally, we will con¬ 
sider exclusively the case of symmetric electrolytes com¬ 
posed of two monovalent ion species, i.e. q + = q_ = 1, 
each species with equal bulk densities p+b = P-b = Pib- 
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III. RESULTS 


A. Solvent configuration 


We investigate herein the partition of solvent molecules 
in the slit pore. To this aim, we will derive the number 
density and the average polarization of solvent molecules 
confined to the pore. We note that the dipolar part of 
the liquid grand-potential was derived in Ref. [5] in the 
form 

= -p sb [ dr ^ e -^(r,a) e -Q[0(r)-0(r+a)] (g) 

J 4tt 

xe -^-5t) d (r,a) 


where the multiple integral is taken over the spatial and 
orientational dipole configurations. Namely, the vec¬ 
tor r indicates the position of the positive charge on 
the dipolar molecule. Then, the infinitesimal vector 
dfl = sinf? d0 dp denotes the solid angle characteriz¬ 
ing the dipolar orientation in the spherical coordinate 
system. The projection of the dipolar orientation vec¬ 
tor onto the z —axis is related to the azimuthal angle as 
a z = a cos 9 (see Fig. [l]). Moreover, the dipolar potential 
in Eq. © is u> s (r,a) = w + (r) + u;_ (r + a) + u> c (r + a/2), 
where the functions w±(r) are wall potentials acting on 
the positive and negative charges of solvent molecules, 
and w c {r + a/2) is coupled to their centre of mass. Re¬ 
defining now the position vector by shifting the latter 
to the centre of the solvent molecule as r' = r + a/2, 
evaluating the density from the functional derivative 
Pd{ r') = 8Qd/5w c ( r'), taking into account the planar 
symmetry of the system, and performing a second trans¬ 
formation a z —> t = z — a./2 on the variable of the in¬ 
tegral over dipole rotations, one gets the dipole number 
density in the form 


Ps{z) 


Psb 


r t2(z) ^ 


_ e -^~Sv a (t,2z-2t) 


'ti(z) 


xe 


Q[0(2z-i)-0(i)] 


(7) 


with the integration boundaries 

t\(z) = max( 0,2 — a/2,2z — d) (8) 

t 2 {z) = min(2^, z + a/2, d). (9) 

The discrete form of Eq. |7]) is reported in Appendix |B] 
We also note that in the present work, the ionic and dipo¬ 
lar potentials 8vi{z) and 5v d {z) in the above equations 
will be calculated within the Debye-Hiickel (DH) approx¬ 
imation. This point is also explained in Appendix [Bj 
We will show that solvent partition in the slit is driven 
by a competition between interfacial polarization and 
surface charge effects. The importance of each effect can 
be accurately described by the average dipole fluctua¬ 
tions m [32]. The latter can be expressed in terms of the 
number density 0 as 

(a 2 z ) 

p m (z) = f- (10) 

a 2 Pd{z) 


(a) 



(b) 



FIG. 2: (Color online ) (a ) Solvent number density |t]| and (b) 
average polarization (101 for the charged solution confined in 
a slit of size d = 1 nm and various surface charges displayed 
in the legend of the bottom plot. The membrane permittivity 
and salt concentration are respectively e m = 1 and pu, = 0.01 
M in both figures. The insets zoom on the interfacial region 
between z = 0 A and 2 A. The horizontal curves in (a) and (b) 
mark respectively the reservoir concentration and the limit of 
freely rotating dipoles p m (z) = 1/3. 


with the statistical average 

r* z(z) Af n 2 

(a?) = / - al (11) 

Jt.iz) a 

Xe Q[0(22-t)-0(t)] 

where a z = 2(z — t). We note that the limit p,(z) = 1/3 
corresponds to the case of freely rotating dipoles exclu¬ 
sively subject to entropy [91 [32]. This situation occurs in 
the bulk reservoir or far from the interfaces where image- 
dipole and surface charge-dipole interactions are absent. 
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Then, in the regime p m (z) <1/3 where repulsive image 
charge forces dominate surface charge effects, the former 
aligns dipoles parallel with the dielectric plane. This re¬ 
gion is expected to be characterized by solvent depletion. 
Finally, for p m {z) >1/3 where surface charge-dipole cou¬ 
pling overcomes image-dipole interactions, the charged 
interface attracts the positive charge and repels the neg¬ 
ative charge of each dipole, resulting in the dipolar align¬ 
ment perpendicular to the interface. In this region, the 
surface charge should induce a local enhancement of the 
solvent density. 

We illustrate in Figs. [2^a) and (b) the number density 
and average polarization of dipoles in a slit of size d = 1 
nm. The model parameters are reported in the caption. 
The membrane permittivity set to e m = 1 is character¬ 
istic of membrane nanopores. First, one notices that for 
intermediate to high surface charges, solvent density ex¬ 
hibits structure formation similar to the one observed in 
MD simulations [34( 135| . Furthermore, the comparison 
of the top and bottom plots shows that the dipole den¬ 
sity and the average polarization are strongly correlated, 
with the position of the minima and maxima coinciding 
almost exactly. First of all, in the top plot, one notes 
that for weakly charged membranes a s = 0.1 e nm -2 , 
solvent molecules are depleted from membrane walls over 
the 3 A thick layer. Then, in the bottom plot and over 
the same region, the low values of the order parameter 
p m {z) <1/3 indicates dipolar alignment parallel to the 
plane. Both the exclusion and parallel alignment effects 
are driven by the rotational penalty for dipoles as well as 
image-dipole interactions induced by the dielectric dis¬ 
continuity at the interface. 

In fig.^b), at the surface charge a s = 0.5 e nm -2 , one 
sees that the average polarization exhibits a first peak 
above the free dipole limit. Thus, over this narrow re¬ 
gion where surface-dipole interactions dominate image- 
dipole forces, solvent molecules show a tendency to align 
perpendicular to the charged wall. In the top plot and 
for the same surface charge, the amplified dipole-surface 
coupling is seen to yield as well a small peak in the den¬ 
sity curve. Then, a further increase of the surface charge 
to the value a s = 1.0 e nm' 2 amplifies the first peak of 
the average polarization and dipole density curves, and it 
also gives rise to a second peak in the polarization. This 
corresponds to the surface charge regime where structure 
formation takes place. However, due to the presence of 
strong image forces in the confined pore, the solvent par¬ 
tition is still characterized by an overall exclusion from 
the pore, that is Pd{z) < p s t,. By setting the surface 
charge to the higher value a s = 2.0 e nm -2 , the pro¬ 
nounced surface charge attraction leads to three peaks in 
the density curve where solvent density exceeds the bulk 
value. The first two peaks are separated by a second sol¬ 
vent exclusion layer. Most importantly, at the position 
of the first peak, the solvent density reaches twice the 
reservoir concentration. We emphasize that such a strong 
interfacial enhancement of solvent number densities has 
been observed in MD simulations for both hydrophobic 


and hydrophilic surfaces (see e.g. Fig.l of Ref. [54]). 

In order to investigate the effect of confinement on sol¬ 
vent partition, we plotted in Fig. [3] the solvent number 
density for various pore sizes ranging between d = 3 A 
and 20 A. One notes that the decrease of the pore size 
that amplifies image-dipole interactions intensifies sol¬ 
vent exclusion from the pore. This effect is particularly 
noticeable below the nanometer pore size (i.e. for d < £ w ) 
where the midpore density is significantly lowered. We 
emphasize that the strength of the dielectric exclusion ef¬ 
fect determining the ionic selectivity of the membrane de¬ 
pends on the pore solvent density since the latter sets the 
intensity of the dielectric discontinuity between the mem¬ 
brane and the nanopore. Thus, the significant solvent 
rejection in narrow pores indicates that solvent-implicit 
nanofiltration models (see e.g. Refs. mm assuming the 
same solvent density in the pore and the reservoir are 
quantitatively inaccurate in predicting the ion rejection 
efficiency of subnanometer pores. In Fig. [3] one also notes 
that the position of the first peak resulting from the com¬ 
petition between image-dipole forces, steric effects, and 
surface charge attraction stays weakly sensitive to the 
pore size. Indeed, multiplying the adimensional position 
of each peak with the corresponding pore size, one finds 
that the former is always located at 0.6 A-0.7 A. 

Finally, in order to evaluate the importance of the di¬ 
electric discontinuity, we calculated the solvent density 
for a dielectrically continuous membrane (i.e. e m = e w 
at the pore size d = 10 A). One observes that in the ab¬ 
sence of image-dipole interactions, the solvent decrement 
layer is replaced by an increment layer with p(z) > p s t for 
z > 0.5 A = a/2. At z = 0.5 A, we notice a sharp peak 
where steric effects resulting from finite solvent size a = 1 
A take over surface charge attraction and decrease the 
density of solvent molecules subject to rotational penalty 
below this point. The closeness of the peak positions for 
e m = 1 and e m = e w suggests that in explicit solvent 
MD simulations, the position of the first density peak is 
mainly determined by the finite solvent size rather than 
surface polarization effects. This point is supported by 
previous MD simulations (see e.g. Figs.1(c) and (d) of 
Ref. 02) where the peak positions located at 2.5 A - 
4.0 A for both hydrophobic and hydrophilic interfaces 
roughly correspond to the molecular size of the water 
TIP models. In the light of these results on solvent par¬ 
tition, we will scrutinize next the partition of the liquid 
charge in the nanoslit. 


B. Partition of the liquid charge 

In order to investigate the charge partition in the 
nanoslit, we plotted in Fig. Qa) the cumulative ion den¬ 
sity 

Pcum,i{ Z ) = [ dz' \p+(z') ~ p-(z')\ (12) 

Jo 
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FIG. 3: (Color online) (a) Solvent dipole density rescaled with 
the bulk density for various pore sizes reported in the leg¬ 
end. The horizontal axis displays the position of the molecule 
rescaled by the pore size. The model parameters are the same 
as in Fig. I 2 I with the surface charge a s = 1.0 e nm -2 . Solid 
curves : e m = 1. Dotted red curve : e m = e w . 


obtained from the solution of the explicit solvent PB 
equation 0 (solid navy curve) and the continuum re¬ 
sult (dashed red curve) from the numerical solution of 
the modified PB equation [10] 

d z e(z)d z (j> c {z ) + ^2 qiPi{z) = -a(z), (13) 


with the sharp dielectric jump function e(z) = 
e w 9(z)d(d — z) + £ m [0(— z) + 9(z — d)]. The modified 
PB formalism ignores the non-uniform solvent parti¬ 
tion in the pore. Fig. ([ 4 ]) (a) shows that moving from 
the interface at z = 0 towards the second interface at 
z = d, the cumulative ion density increases monotoni- 
cally and reaches the surface charge value in the mid¬ 
pore p C um,i(d/ 2) = a s and twice the surface charge at 
the second interface, i.e. p C um,i{d) = 2 a s . This is a con¬ 
sequence of the electroneutrality condition according to 
which there must be as many monopole charges in the 
nanoslit as the fixed surface charge on the membrane 
walls. One also notes that the cumulative ion density 
from the NLPB formalism (navy curve) and continuum 
theory (red curve) are very close to each other. Thus, 
solvent structure plays a minor role in the cumulative 
ion charge configuration. 

We now note that integrating Eq. © from the inter¬ 
face at z = 0 to any point z located in the pore, one can 
express the continuum electric field E c (z) = d z (f> c {z) as 


E c {z) = 47r£ w [<j s - Pcum,i{z)\ ■ 


(14) 


We reported the electric field profile E c (z) in Fig. 0b). 
As suggested by Eq. (14) and in agreement with the ion 



FIG. 4: (Color online) (a) Cumulative charge density of sol¬ 
vent molecules p C um,s{z) (solid blue curve), ions pcum,i(z) 
(solid navy curve), and total cumulative charge density 
pc.um,s{z) + p C um,i{z) (dashed-dotted black curve). The 
dashed red curve displays the cumulative ion density of the 
continuum theory, (b) Electric field from the explicit sol¬ 
vent (solid navy curve) and continuum theory (dashed red 
curve). The inset displays the local solvent charge density 
p s +(z ) — p s -(z). In (a) and (b), membrane permittivity is 
£ m = 1, bulk ion density pu, = ICC 2 M, membrane surface 
charge <7 S = 1.0 e nm -2 , and pore size d = 2 nm. 


density curve in Fig. 4|a), the electric field evolves from 
E c (0) = 4n£ w a s at the interface to zero in the midpore 
where p cum ,i(d/ 2) = cr s and the field reaches the value 
E c (d) = —4Tr£ w <j s at the second interface where the cu¬ 
mulative ion density is twice as large as the wall charge. 
In other words, within the continuum formalism, the spa¬ 
tial variation of the field is solely determined by the ion 
screening effect. 
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We consider next the cumulative solvent density 

Pcum, s (z ) = / d z' \p s +{z') - p s -(z')] 

Jo 


(15) 


displayed in Fig. 0a) by the solid blue curve together 
with the local solvent density shown in the inset of 
Fig. 0b). In the latter figure, one notes that the interfa¬ 
cial regions are each characterized by a positively charged 
solvent layer separated by a negative solvent layer. In 
Fi g- 0a), this local charge distribution is shown to result 
in an initial increase of the cumulative solvent density up 
to the maximum value p C um,s(z ) ~ a s where it starts to 
decrease and vanishes in the mid-pore. For z > d/ 2, the 
cumulative density becomes negative and drops mono- 
tonically up to the vicinity of the second interface located 
at z = d where it reaches its minimum value. Moving 
closer to the second interface, due to the interfacial posi¬ 
tive solvent charge layer, the cumulative density rises and 
becomes zero on the pore wall. The latter point shows 
that solvent molecules bring no contribution to the total 
charge in the slit. 

Before considering the total cumulative charge density, 
one notes that by integrating Eq. Q, one gets the equiva¬ 
lent of Eq. (14) for the electric field in the solvent-explicit 
liquid, 

E(z) = 4 iris [<r a - Pcum,i(z) - Pcum,s{z)\ . (16) 


Eq. (16) indicates that the electric field profile is deter¬ 


mined by the combined effects of salt screening (the sec¬ 
ond term on the right-hand-side) and dielectric screening 
associated with polarization charges (the third term). We 
plotted in Fig. 0a) the total charge density pcumXz) + 
Pcum,s(z) (dotted-dashed curve) and reported the elec¬ 
tric field in Fig. 0b) (solid navy curve). Observ¬ 
ing both figures, one first notes that the interfacial sol¬ 
vent depletion responsible for a dielectric screening defi¬ 
ciency results in a surface field significantly larger than 
the field of the continuum theory. Then, due to the rise 
of the cumulative solvent density, the total cumulative 
charge density increases at the interface. As a result of 
the amplified dielectric screening effect, the electric field 
decreases rapidly from the surface value E( 0) = 47 t£b& s 
to the magnitude of the continuum field E c (z). More¬ 
over, we see that far away from the interfaces, the total 
cumulative charge density stays very close to the wall 
charge. Finally, we note that the field E(z) exhibits fluc¬ 
tuations around the continuum field. This arises from 
the solvent charge density fluctuations displayed in the 
inset of Fig. 0b). In the next part, we will investigate 
the effect of non-uniform solvent charge distribution on 
the partition of ions in the nanoslit. 


C. Solvent effects on counterion partition 

In this part, we probe solvent structure effects on the 
partition of counterions in the slit pore. We illustrate 



z(A) 



z(A) 

FIG. 5: (Color online) (a) Rescaled counterion densities in 
the slit of size d = 2 nrn and different charges increasing 
from bottom to top : cr s = 0.7 e nm -2 (black curve), 1.0 e 
nm -2 (navy curve), and 1.3 e nrn -2 (red curve). The inset 
zooms on the interfacial region. We also display the coun¬ 
terion density from the solvent-implicit theory (dashed navy 
curve), (b) Image charge potential Svi(z) (blue curve), elec¬ 
trostatic potentials <f> c (z) from the implicit solvent formalism 
(black curve) and (j>(z ) of the explicit solvent approach (red 
curve) at the pore charge <r s = 1.0 e nm -2 . Reservoir ion den¬ 
sity and membrane permittivity are pib = 0.01 M and e m = 1 
in both figures. 


in Fig. 0a) cation densities from the solvent-explicit for¬ 
mulation (solid curves) for increasing pore charge from 
bottom to top. We also reported the cation density pro¬ 
file from the continuum formulation of Eq. ( 13 ) (dashed 
curve at the surface charge a s = 1.0 e nm -2 ). Fig. 0b) 
displays in turn the electrostatic potentials <p(z) and 
4> c (z) as well as the ionic self-energy Si’i(z)/ 2 and the to¬ 
tal ionic potential of mean force (PMF) (f>(z) + dvi(z)/ 2 . 
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We also note that in this figure, the strong amplitude of 
the repulsive ionic self-energy is known to result from im¬ 
age charge interactions between ions and the membrane 
of low dielectric permittivity. 

In Fig. [5^a) , the result from the continuum theory at 
the surface charge a s = 1.0 e nm^ 2 (dashed navy curve) 
shows that due to the competition between the repul¬ 
sive image-charge potential and the attractive electro¬ 
static potential, the density curve exhibits at both in¬ 
terfaces a counterion depletion over 1 A followed by a 
maximum. The maxima are separated in turn by a well 
in the mid-pore area. At the same surface charge, the 
result of the solvent-explicit theory exhibits the same 
counterion peak, with the corresponding curve exhibit¬ 
ing fluctuations around the continuum result (see the in¬ 
set). Most importantly, at the surface charge a s = 0.7 
e nm~ , another counterion peak located at a closer dis¬ 
tance from the interfaces appears. Increasing the surface 
charge from this value to a s = 1.3 e nm -2 , the height of 
this peak rises and exceeds that of the second peak. The 
additional interfacial counterion adsorption layer absent 
in the continuum theory is a consequence of the non- 
uniform solvent configuration. Indeed, Fig. [5jb) shows 
that due to the interfacial dielectric screening deficiency 
discussed in Section |IIIB[ the potential of the solvent- 
explicit formalism <f>(z) is much larger than the contin¬ 
uum potential (f> c (z ) in the interfacial regions. More pre¬ 
cisely, at the surface charge a 3 = 1.0 e nm -2 , the surface 
potential values are 0 C (O) ~ —8 and 0(0) ss —27 
ksT. In other words, close to the pore walls, the am¬ 
plitude of the surface charge induced electrostatic force 
becomes comparable with image-charge forces. The lo¬ 
cation of the adsorption peaks absent in the continuum 
approach corresponds precisely to the region where the 
repulsive image-charge force is locally dominated by the 
strongly attractive electric field. We emphasize that the 
presence of this interfacial counterion adsorption layer 
may have important consequences in the functioning of 
nanofluidic devices [34] [35]. We plan to explore this issue 
in an upcoming article. 

We scrutinize now the effect of the membrane polarity 
on the ionic structure formation illustrated in Fig. Ha). 
To this aim, we plotted in Fig. [6] counterion densities at 
three different membrane permittivities. The main plot 
shows that increasing the membrane permittivity from 
£ m = 1 to 30, the reduction of image charge forces am¬ 
plifies the height of the first counterion peak. Because the 
electroneutrality condition fixes the total number of ions 
in the pore, this amplified counterion adsorption is com¬ 
pensated by a density decrease in the mid-pore. More¬ 
over, with a further increase of the membrane polarity 
to the value e m = e w where image-charge forces vanish, 
the two counterion peaks merge and leave a shouldering 
located at a distance ~ 1 A from each interface. In this 
limit where the ionic structure disappears, the counterion 
density drops monotonically from the interface towards 
the mid-pore. We note that this behaviour is characteris¬ 
tic of the MF counterion densities previously investigated 



z(A) 


FIG. 6: (Color online) Counterion densities in the slit of size 
d = 2 nm and surface charge a 3 = 0.7 e nm' 2 for various 
membrane permittivities : e m = 1 (black curve), e m = 30 
(red curve), and e m = e w (blue curve). The inset shows the 
net electrostatic force F to t(z) = —d z [4>(z) + 5vi(z)/ 2] acting 
on each counterion. 

in Ref. [24] for single charged planes. 

Finally, in terms of the change in the membrane per¬ 
mittivity, we will show that the ionic structure for¬ 
mation results from the competition between surface 
charge attraction and image-charge repulsion. To this 
aim, we plotted in the inset of Fig. [6] the electrostatic 
force F to t{z) = —d z \(j){z) + Svi(z)/ 2] experienced by each 
counterion located in the nanoslit. One sees that at the 
membrane permittivity £ m = 30 (red curve), the net force 
switches several times from repulsive ( F tot {z) > 0) to at¬ 
tractive ( F to t{z ) < 0). At the points where the net force 
vanishes ( F tot {z ) = 0), surface charge attraction com¬ 
pensates exactly image charge repulsion. This gives rise 
to the counterion peaks displayed in the main plot. In 
the limit £ m = £ w where structure formation disappears 
(blue curve), the net force is seen to be purely attractive 
in the nanoslit ( F tot {z ) < 0). This leads in turn to a 
monotonical counterion increase towards the wall. 


IV. CONCLUSIONS 

In this work, we investigated the coupled effects 
of membrane charge, polarization forces, and non- 
uniform dielectric permittivity on the partition of solvent 
molecules and ions under nanoconfinement. We showed 
that the solvent configuration is characterized by a strong 
correlation between dipolar orientation and solvent num¬ 
ber density. Namely, local solvent excess is usually ac¬ 
companied with the alignment of solvent molecules per¬ 
pendicular with the pore walls, whereas solvent exclu¬ 
sion layers are associated with dipolar alignment parallel 
















with the membrane. We also found that even in strongly 
charged membranes, subnanometer pores are character¬ 
ized by a reduced solvent density and dielectric permit¬ 
tivity. This indicates that solvent-implicit ion rejection 
models assuming the same pore and reservoir permittiv¬ 
ities are not reliable in predicting the ionic selectivity of 
subnanometer pores Hill- Moreover, the non-uniform 
polarization field of solvent molecules resulting in an in¬ 
terfacial dielectric screening deficiency leads to a strong 
surface electric field. This results in a counterion adsorp¬ 
tion peak absent in continuum theories m- The obser¬ 
vation of this additional counterion layer associated with 
solvent structure is the main result of our work. This pe¬ 
culiarity may have important effects on nanofluidic ion 
transport [Ml[55] , 

The present work aimed at probing the equilibrium 
properties of confined electrolytes associated with sol¬ 
vent charge structure. In an upcoming work, we plan 
to explore the impact of the points investigated herein 
on nanofluidic transport. It is also important to point 
out the limitations of the present approach. First, we 
emphasize that our non-MF NLPB formalism is hybrid 
since we approximated the electrostatic propagator ac¬ 
counting for charge correlation effects with the solution 
of the continuum DH equation. The reason for this ap¬ 
proximation is the absence of any recipies for the solution 
of the solvent-explicit non-local kernel equation derived 
in Ref. f(5]. Moreover, although our formalism accounts 
for the extended solvent charge structure, the associated 
scalar field theory is a continuum formulation since it 
does not include excluded volume effects and dielectric 
cavities associated with solvent molecules. Despite these 
limitations, our approach goes beyond the classical con¬ 
tinuum theories of electrolytes since it is the only formal¬ 
ism that can qualitatively reproduce interfacial ion and 
solvent structure formation observed in MD simulations. 
The comparison of the emerging ion transport properties 
with experimental data is of course needed in order to 
check the accuracy of the present theory and the conse¬ 
quences of the above-mentioned limitations. 


and the Bessel function of the first kind Jq{x) [Mi- In 
Ref [9], we had derived a variational kernel equation for 
the electrostatic Green’s function u(r, r'). At present, 
an analytical or numerical approach to solve this highly 
non-local and non-linear kernel equation is not available. 
Thus, we will approximate the electrostatic Green’s func¬ 
tion in Eqs. (A1 )-( [A2| ) by the solution of the local DH 
equation in slit pores [2J, 


v(z,z',k) = 


2ttL, 


{ C 


-p\z-z'\ 


(AS) 


A 


g -p(z+z') _|_ e p{z+z'-2d) 


1 - A 2 e~ 2 P d V 

+2Ae" 2pd cosh (p\z — /()]} 


where we have introduced the Bjerrum length in water 
solvent t w = e 2 /(Ane w kT) and the functions 


P = 
A = 


a/k 2 + k 2 
^wP ^rr;A 

£wP T ^rnk 


(A6) 

(A7) 


with the DH screening pa rame t er n 2 = 4 ir£ w QiPib- 
By substituting into Eqs. (Al)-(A2) the Fourier trans¬ 
formed propagator (A5), the explicit form of the ionic 
and dipolar self-energies read 

d kk A 


Svi(z) = t w f 
Jo 

5v d {z,a z ) = £ w 

Jo 


p 1 - A 2 e -2pd ^ 

x je _2pz + e - 2p{d ~ z) + 2Ae _2pd j 
3 d kk A 

7wAVV f(z ' az) ’ (A9) 


where we introduced the auxiliary function 
F(z,a z ) = e~ 2pz + e - 2 P( d - z ) + e -Md-*-a.) ( Al0 ) 

+e -2 P (z+oj + 4Ae -2 pd 

-2 Jo [fc|o,| I] j e - p ( 2 *+“*) + e —J>(2d-2*-o.) 

+2Ae _2pd cosh(pa z )} . 


Appendix A: Derivation of the ionic and dipolar 
self-energies 


Appendix B: Relaxation algorithm for the solution 
of the non-MF NLPB equation ([l]) in slit pores 


We report in this appendix the ionic and dipolar self¬ 
energies in Eqs. § and 

Svi(z) = v(z, z) — v b (0 ) (Al) 

Sv d (z, a z ) = v d (z, a z ) - 2e b (0) + 2 v b (a), (A2) 


with the electrostatic Green’s function and dipolar po¬ 
tential 

r°° <\kk r i 

u(r,r') = j( — J 0 [fc|r || -rf | |]u(z,A (A3) 

r°° 

Vd(z,a z ) = J — {v(z,z) + v(z + a z ,z + a z ){M) 
-2 v(z,z + a z ) J 0 [fc|a|||] } , 


We present herein a numerical relaxation scheme for 
the solution of Eq. (JT|) in slit pores. The scheme will be 
an extension of the algorithm presented in Ref. 9] for 
simple interfaces to the case of slit pores. In order to 
simplify the form of Eq. ([l]), we note that we considered 
in the present work the case of solvent molecules with 
monovalent elementary charges Q = 1 and symmetric 
electrolytes composed of two monovalent ion species, i.e. 
q+ = q- = 1, with the bulk ion densities = p-b = Pib- 
With these simplifications, the NLPB Eq. (|T|) takes the 
form 

d 2 z ^{z) - £ W K 2 e ~ 5vi ^ /2 sinh [<j>(z)] 


(Bl) 
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(•a 2 ( 2 ) 


'ai(z) 


d a z 
2 a 


e -5v d (z,a,)/2 ginh ^ + az) j 


= 47T £ B a s [ 6 ( 2 ) + S(d - z )\, 

where we introduced the solvent screening parameter 


k 2 = 8n£ B p s b- By integrating Eq. (Bll in the vicinity 


of the interfaces at 2 = 0 and 2 = d, one obtains the 
two boundary conditions required for the solution of this 
equation, 


A'(0) = 4rf B (T s 
f>'(d) = -Ant B cra- 


(B2) 

(B3) 


The relaxation algorithm is based on the discretization 
of Eq. ( |B2[ ) on the 2 axis. This is done on a discrete lattice 
located between 2 = 0 and 2 = d. The lattice is composed 
of 2 N + 1 mesh points with separation distance e. In 
order to express Eq. (Bll on the lattice, we introduce 
the discretized coordinate z n = e(n — 1) and potential 
4>n = (f>( z n) with 2N + 1 < n < 1 . By introducing as well 
the discrete form of the Laplacian operator e 2 (j>''(z) = 
ip n + 1 + V’n-i — 20„, Eq. (Bl) can be expressed on the 
lattice as 

■0 n = \ {0n+i + 0n-i - r e~ 5vi(Zn)/2 smh(i/j n ) (B4) 
h (n) 

-s e- 5vd ^’ a ^- n + l)/2 smh (0„-0j) 

j=ji(n) 

with the parameters r = E w e 2 K 2 , s = e w £ 3 k 2 , and the 
functions 

ji(n ) = max(l, n — n a + 1) (B5) 

j 2 (n) = min(2AT + 1, n + n a — 1), (B6) 


where the index n a is defined as z(n a ) = a . Eq . (B4) will 
be solved with the boundary conditions (B 2 )-(B3) that 


can be expressed on the lattice as ipo = 0 i — 4'k£ b <J s £ 
and 02AT+2 = 02iv+i ~ 4ir£ B cr s E. We finally note that on 
the same discrete lattice, the solvent number density 0 
takes the form 


h(r, 


Pd(z n ) = - > 

a z— 1 

j=h (n) 


g— ^Svd(Zj ,a z ^2z„-2zj) — 


(B7) 


with the boundaries 


h(n) = max(l, n — n' a + 1,2n — 2iV — 1) (B 8 ) 

I 2 (n) = min(2n — 1, n + n' a — 1,2JV + 1), (B9) 

and the index n' a defined as z(n ' a ) = a/ 2 . 


The relaxation algorithm consists in solving Eq. (B4) 


by iterating an initial reference potential (/> r {z). As al¬ 
ready explained in Ref. [9], the complication stems from 
the fact that for the convergence to be achieved, the in¬ 
put function has to verify the same boundary conditions 


as the NLPB Eq. (Bl), i.e. the guess potential should 
obey Eqs. (B2 |-( |B3[ )~ These conditions are not verified 
by the local MF PB equation since the PB potential sat¬ 
isfies a different boundary condition 0'(O) = 47T £ w a s at 
2 = 0. In other words, the latter does not take into 
account the dielectric void in the neighbourhood of the 
charged interfaces. 

In order to account for the interfacial dielectric screen¬ 
ing deficiency absent in the PB formalism, we will reit¬ 
erate the trick explained in Ref. [9] and extend the guess 
potential derived in this earlier work for simple interfaces 
to slit nanopores. Namely, neglecting the dipolar self¬ 
energy dv(z, a z ) in Eq. (Bl) and linearizing this equation 
in terms of the potential 4>{z), one finds that the result¬ 
ing equation is reduced for 2 < o and e w k 2 <C k 2 to 
4>"{z) — c 2 K 2 (j){z) = —47r£scr(2), where we introduced the 
geometric factor c 2 = 1/2 associated with the rotational 
penalty for dipoles in the vicinity of the rigid interfaces. 
The solution of this equation yields the electrostatic field 
in the form </>' (z) = 47r^Btr s e _CKs2: . We now note that 
this field is also the solution of the non-uniform Pois¬ 
son equation d z e{z)d z (j>{z ) = — An£Bcr(z), with the di¬ 
electric permittivity function given by e(z) = e CK,sZ . In 
Ref. 0 , this exponential law was shown to reproduce ac¬ 
curately the behaviour of the dielectric permittivity up 
to the characteristic distance g?i = ln(e UJ )/(cK s ) « 0.3 A 
where the non-local permittivity reaches the bulk value 
(see Fig.6(a) of Ref. [5])- Inspired by this observatio n, we 
approximate for pores with thickness d > 2d\ Eq. (Bl) 
by the following equation, 

d z e{z)d z <j>{z) - £ w k 2 c {z) [c, b r {z) - 0 O ] = -47r£ B er(2), 

(BiO) 

where 0o is a uniform Donnan potential that will be 
determined from a variational minimization procedure. 
Furthermore, the piecewise ionic screening parameter 
reads k c (z) = k9(z — d\)0{d — c?i — 2 ), and the inho¬ 
mogeneous dielectric permittivity function is given by 

e(z) = e CK,sZ 9{di — z)d(z) + e w 0{z — d\)9(d — d\ — z) 

e CK °(d-z) 9 ( z _ d + dl } 9 ( d _ z y ( B11 ) 


The solution of Eq. ( |B10 ) satisfying the continuity of 
the potential 4>(z) and the displacement field e(z)(j)'(z) 
at 2 = d\ and 2 = d — d\ reads 


Mz) — 00 


A 


47t£ B (Ts c 


CK S 


9{di - z)9(z) (B12) 


47r£ w a s cosh [re(d /2 — 2 )] 
k sinh [n{d/2 — di)] 

A + 4 rf B^ c - CK ,(d- Z ) 


CK, 


where we introduced the constant 


9{z — di)6(d — di — 2 ) 


9{z — d + di)6(d — 2 ), 


A = 47r£,„cr, 


coth 


K I ^ ^ d i 


1 

CK S 


(B13) 


The solution scheme consists in injecting into Eq. (B4) 


the guess potential (B12) at the first iterative step, using 
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the output function as the new input potential at the next 
iterative step, and continuing the cycle until numerical 
convergence is obtained. 

Finally, we explain the variational determination of 
the Donnan potential <f >o that takes into account non- 
linearities neglected by the linear trial form (B12). The 


By taking the derivative of Eq. (B14) with respect to the 


Donnan potential cj> o, one obtains 

rd 


— 2 pu, f d ze Svi ( z )/ 2 sinh[^ r (^)] — 2cr s = 0. (B15) 

Jo 


approach consists in injecting the guess potential (B12) The integral equation (B15) can be easily solved with 


into the part of the variational Grand potential that gen¬ 
erates Eq. ([!]). Up to constant factors that do not depend 
on the variational parameter (f> o, the variational func¬ 
tional to be optimized with respect to (j> o reads [5] 

fd 

h[<po\ = ~2pib / dz e~ Sv d z )! 2 cosh[0 r (jz)] - 2 a s <j) 0 . 

Jo 

(B14) 


a standard dichotomy algorithm in order to obtain the 
numerical value of the Donnan potential. 
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